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Abstract 

In this paper, we investigate a Langevin model subjected to stochastic intensity noise (SIN), 
which incorporates temporal fluctuations in noise-intensity. We derive a higher-order Fokker- 
Planck equation (HFPE) of the system, taking into account the effect of SIN by the adiabatic 
elimination technique. Stationary distributions of the HFPE are calculated by using the perturba- 
tion expansion. We investigate the effect of SIN in three cases: (a) parabolic and quartic bistable 
potentials with additive noise, (b) a quartic potential with multiplicative noise, and (c) a stochastic 
gene expression model. We find that the existence of noise intensity fluctuations induces an intrigu- 
ing phenomenon of a bimodal-to-trimodal transition in probability distributions. These results are 
validated with Monte Carlo simulations. 
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1. INTRODUCTION 



Many real- world systems are inhomogeneous and fluctuant. Stochastic processes are 
often used for modeling such fluctuant systems in many fields, including physics, biology, 
and chemistry. The dynamics in these systems can be described by a Langevin equation 
given by 

(IT 

- = f(x) + g(x)Ut), (1) 

where f(x) = —d x U (x), U(x) denotes a potential, g(x) is an arbitrary function of x, and £ x (t) 
is the white Gaussian noise with correlation {£, x {t)!; x (t')) = 2D x S(t—t'). Although white noise 
can reflect microscale properties of fluctuations, it is uniform when seen from mesoscopic or 
macroscopic time scales [Fig. Ufa)]. One widely used approach for incorporating mesoscopic 
or macroscopic inhomogeneity in noise is colored noise, where the white noise £ x (t) in Eq. (TTJ) 
is replaced by z(t) with the Ornstein-Uhlenbeck process: 

dz z £ z (t) ,„ s 
at t t 

where £ z {t) expresses the white Gaussian noise [(£z(t)£z(t')) = 2D z 5{t — t')). Equa- 
tion (j2J) yields colored noise with the finite correlation time r, that is, (z(t)z(t')) = 



(D z /t) exp { — \t — t'\/r}. The existence of corre 



ation time in the noise sources can induce 



many phenomena, such as resonant activation [1] and noise-enhanced stability {2 — 6 1 . 

In the present paper, we deal with mesoscopic time-scale inhomogeneity in a way other 
than with colored noise; here, we consider temporal noise-intensity fluctuations. We assume 
that the noise intensity in Eq. (CQ) is not constant but modulated by the Ornstein-Uhlenbeck 
process. Our model is described by the following coupled Langevin equations instead of 
Eq. (JO: 

n T 

^ = f(x)+g(x)sUt), (3) 
ds 

_ = _ 7 ( s _ a ) + v% (f). (4) 

Here, £ g (t) is the white Gaussian noise [(£ s (0£s(O) = 2D s 5(t — t')], and 7 and a are the 
relaxation rate and the mean of the Ornstein-Uhlenbeck process, respectively. The intensity- 
modulated noise term s£ x (t) in Eq. ([3]) is herein called stochastic intensity noise (SIN) [Fig. 
1(b)]. This point of view was originally introduced in Heston's stochastic volatility models 
in financial engineering |7| and has since been analyzed in econophysics [8]. With f(x) oc x 



and g(x) oc x, Eqs. fl3]) and (SJ) are similar to those in the Heston model, where the variance 



he 



9] 



of noise is governed by the Feller process (also referred to as the square-root process or t 
Cox-Ingersoll-Ross process). Escape events in the Heston model were studied in Ref. 
using a cubic potential. Note that the variable s in Eqs. (EJ) and (jl]) takes negative as 
well as positive values since our model can be considered as white noise with multiplicative 
term s, which is in contrast with the positive variance in the Heston model [7]. In physics, 
superstatistics includes the concept of temporal and/or spatial environmental fluctuations 
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1^ | . Superstatistics was originally introduced, under specific conditions, to account for 
asymptotic power-law distributions (e.g., g-exponential distributions and q- Gaussian distri- 
butions) that emerge as maximizers of non-additive (Tsallis) entropy 15j,ll6[|. Superstatistics 
has since been applied to the interpretation of quasi-equilibrium thermodynamics, and con- 



cepts of superstatistics have also been applied to stochastic processes [111 . 1 131 . |17H20|. In 
particular, Ref. [l7J extended superstatistics to the path-integral representation and showed 
that some stochastic models can be covered by this representation. A direct connection 
between Tsallis statistics and financial stochastic processes was indicated in Ref. [21 1. 

In many biological and chemical processes, the relaxation time of x may be larger than 
that of environmental fluctuations (noise-intensity processes). This is the case in stochastic 
gene expression models in which the decay time of x is on the order of minutes 22J (Sec. 15 .3) 1 . 
Bearing this fact in mind, we will investigate systems driven by SIN with faster decay time 
compared with that of x (7 ^> 1). Furthermore, in real- world systems, we expect that the 
f(x) in Eq. ([3]) is often given by a complex nonlinear function and is also accompanied 
by nontrivial multiplicative noise expressed by an appropriate g(x). These properties are 
different from those for financial engineering. In order to obtain the probability distribution 
P(x, t), we use adiabatic elimination with eigenfunction expansion 23J. We obtain a higher- 
order Fokker-Planck equation (HFPE) with higher-order derivatives, not included in the 
conventional Fokker-Planck equation (FPE). We calculate stationary distributions of the 
systems described by the HFPE by using the perturbation expansion. 

To investigate the effects of SIN, we consider stationary distributions for three cases: 
(a) parabolic and quartic bistable potentials with additive noise, (b) a quartic potential 
with multiplicative noise and (c) a stochastic gene expression model. In the additive noise 
case (a), we show that SIN changes the stationary distribution from exponential forms (the 
Boltzmann-Gibbs distribution) to non-exponential forms. At the same time, the stationary 
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FIG. 1: (Color online) Intuitive paths of (a) white noise and (b) SIN whose intensity is governed 
by the Ornstein-Uhlenbeck process [Eq. Q]. 

distributions of case (a) are sharpened because of noise- intensity fluctuations. In case (b), we 
show that the existence of noise-intensity fluctuations induces the transition of distributions: 
the stationary distributions are uni-, bi-, or trimodal, depending on the SIN parameters. It 
is important to note that the trimodal distribution does not emerge under white noise. In 
case (c) which we call the gene expression model, a nonlinear function is given as a drift 
term (change in expression levels). Thus, case (c) shows that our approximation scheme can 
be applied to general configurations that include non-trivial drift terms. 

The remainder of this paper is organized as follows: In Sec. |2} we describe the model 
proposed in this paper. Adiabatic elimination with eigenfunction expansion is applied to 
the model in Sec. [3] (details of the derivation of the HFPE are explained in the Appendix) . In 
Sec.Hl we proceed to the calculation of the stationary distributions of the obtained HFPE by 
using of the perturbation expansion. In Sec. El we investigate effects of SIN in the three cases 
(a), (b), and (c) mentioned above. In Sec. El we analyze effects of higher-order derivatives 
in the HFPE on positivity and moments of distribution functions, and also consider the 
opposite case, in which a decay time of the noise-intensity fluctuations is very slow (7 — > 0). 
Finally, we give the conclusions in Sec. [3 

2. THE MODEL 

We consider the Langevin equations given by Eqs. fl3]) and (JIJ. Since Eq. (@J does not 
depend on x, the stationary distribution P st (s) of s is easily obtained: 
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Time evolution of the probability distribution P(x, s; t) is given by 

d 

—P(x,s;t) = L FP (x,s)P(x,s;t), (6) 
where Lpp(x, s) is a FPE operator composed of 

L FP (x,s) = L x (x,s) + 7L s (s), (7) 

with 

L,(x, s) = -— {f{x) + D x g{x)g'{x)s 2 } + D x s 2 —g{x)\ (8) 

L.(s) = ±(8-a) + D.£ 5 . (9) 

We employ Stratonovich's calculus because it is expected to be more relevant to physical 
applications than Ito's [241 ] . We calculate a projected time evolution equation of P(x;t) = 
f ds P(x, s;t) that can be applied to general configurations. 



3. ELIMINATION OF THE FAST VARIABLE 



In this section, we eliminate s from Eq. ([6]) by using adiabatic elimination under the 
assumption that 7 ^> 1 (i.e., the decay time of s is much faster than that of x). Adia- 
batic elimination has its origin in the Born-Oppenheimer approximation. Because of its 



usefulness, adiabatic elimination has been applied to many stochastic systems 
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25 
28 



26j. In 



These 



stochastic volatility models, 7 — > adiabatic elimination has been applied 
models take advantage of the fact that noise intensity (volatility) changes on a macroscopic 
time scale in financial markets (this description agrees with Beck's superstatistical Brown- 
ian motion [13]). In these studies, the obtained solutions do not explicitly include 7 as a 
parameter. In the present paper, we consider the 7 ^> 1 case and derive a HFPE containing 
0(7 _1 ) terms by using adiabatic elimination with eigenfunction expansion 23]. 

By applying adiabatic elimination up to 0(7 _1 ), we obtain the following HFPE (see the 
Appendix for the derivation): 



d_ 

Of 



P(x;t) 



Af( x ) + QA g + ^A 2 g }P(x;t), 



with 



A, 



dx 2 



g{x) 



d_ 

dx 



g'(x)g(x), 



(10) 



(11) 



where Q and R in Eq. ( fit)]) are defined as follows: 

Q = D x (D s + a 2 ) , 

R = D\D S (4a 2 + D s ) . 
For the case of additive noise [i.e., g{x) = 1], Eq. flTDJ) reduces to 

|p(, ;i) = {-| /w+Q ^ + ^}p(, ;t ). (12) 

Because Eqs. (|T0|) and (fl2|) have derivatives of orders higher than two, Eqs. ({TO]) and (|T2|) 
are referred to as the HFPE. In systems driven by colored noise, adiabatic elimination up 
to O(t) can be expressed by the conventional FPE [23J [see Eq. (J2J)]. On the other hand, 
0(7 _1 ) terms are accompanied by non-FPE terms in our SIN case. Therefore, in order to 
incorporate the effect of 7, we must use the HFPE form of Eq. fflUj) . For 7" 1 — > 0, the last 
term of Eq. ( TTOj) vanishes, and we obtain FPE: 

-P(x; t) = -— {/(*) + Q</(xM*)} P(x; t) + Q—g(x) 2 P(x; t). (13) 

Eq. f)13p corresponds to conventional (7 — > 00) adiabatic elimination. Q plays the role of 
the effective noise intensity of the corresponding white Gaussian noise process. 

4. STATIONARY DISTRIBUTION 

In many practical cases, stationary distributions play important roles. In this section, we 
calculate the stationary distribution P s t{x) of Eq. (fTOj) . which yields the following differential 
equation: 

/ r, ^ , ( 9 „2 j \ D ,s , / ^ o , \ d ( 8 2 



t Pst(x) + " ^ J Pst(x) + £ I ^ - 9,9 ) Tx W - 9,9 > P *M = °' (14) 



= J R D a £ 8 (4a 2 + £g 
5 7 Q 7(^ s + « 2 ) ' 1 j 



where 



It is easy to see that a solution of Eq. ( 1141) for e = is the stationary distribution of 
Eq. ( fl3l) . Thus, for e <C 1, it is expected that a solution of Eq. ( Tl4"l) can be approximated by 
the perturbation from the stationary distribution for e — 0. In solving Eq. (1141) . we adopt 



the perturbation expansion 29l432| . given by 



P st {x) = n (x) + eU^x) + e 2 U 2 (x) + ■ ■ 



We specifically define the following truncated first-order approximation: 

P$l\x) = U (x) +sU 1 (x), (16) 

where n (x) is the stationary distribution of the unperturbed case [Eq. (ITS]) ], and H\(x) 
corresponds to the first-order correction term. Substituting Eq. ( JT6"j) into Eq. (|14p and 
comparing the order of e up to 0(e), we obtain 

0(1) ! [ ^- + g'(x)g(x)^n (x) = ^g(xrn (x), (17) 

0(e) |M + ^ (x) |n 1 (a ; ) = ^(x) 2 n 1 (a;) + 0(x), (18) 



where 



<>(.r) <{ -^-g(x) 2 - g'(x)g(x)\ ^- <-^g(x) 2 -g'(x)g(x) } Ll u (r) 



From Eq. (IT7|) . n (x) is given by 



^o( x ) = e x P ( -pz I dv 



Z\g(x)\ ^\QJ g(v)\ 

where Z is a normalizing term [J rfxllo(x) = 1]. Since Eq. ( JT8l) is a first-order differential 
equation, Ili(x) can be obtained analytically in many cases, as will be discussed in Sec. El 
The perturbation expansion yields reliable results for e < 1. Combining the approxima- 
tion condition of adiabatic elimination and the perturbation expansion, it can be concluded 
that we are able to calculate the stationary distribution by Eq. (fl6]) for systems with 7> 1 
and e <C 1. 



5. MODEL CALCULATIONS 



We apply the approximation method obtained in Sec. H] to three cases: parabolic and 
quartic bistable potentials with additive noise (Sec. 15 . 11) . a quartic potential with multi- 
plicative noise (Sec. 15. 2p . and a stochastic gene expression model (Sec. 15.31) . We analyze 
effects of noise-intensity fluctuations on shapes of the stationary distributions. We also car- 
ried out Monte Carlo (MC) simulations to determine the reliability of our approximation. 
In the following, we show results of three calculation methods: 

• HFPE stationary distribution 

The stationary distribution of the HFPE [Eq. ffTUl) ] is shown using the perturbation 

7 



expansion explained in previous sections including the correction term Ui(x). This is 
an 0(7" 1 ) approximation. 

• FPE stationary distribution 

The stationary distribution of the FPE [Eq. (Tl3~j) ] is shown. This is identical to Ilo(x) 
and does not include the effect of noise-intensity fluctuations. This case corresponds 
to 7~ 1 —7-0 adiabatic elimination. 

• MC simulation 

We have employed a simple Euler-forward scheme with a time s tep size At = 10~ 6 
(for details of the algorithm, readers may refer to Sec. 3.6 of Ref. 33J). We have used 
A = 5 x 10 7 points to calculate empirical distributions. For all the line-symmetric 
potentials, MC data are symmetrized with respect to x = 0. 

For evaluating quantitatively the reliability of HFPE stationary distributions, we calculate 
the root-mean-square (RMS) distance between the HFPE (or FPE) and MC distributions. 
We first divide an interval of x (window of each figure) into M points (M = 100), each of 
which is denoted by Xi. Let Pi and P s t(xi) be density values of the MC and HFPE (FPE) 
distributions at Xi, respectively. RMS is defined by 



RMS 



M 

^Mg( P "- F " M )' 



5.1. Parabolic and Quartic Potentials with Additive Noise 

We investigate effects of additive SIN in parabolic and quartic bistable potentials. We 

first calculate the stationary distribution of the parabolic potential [U(x) = x 2 /2] with 

additive SIN [g(x) = 1]. According to Eqs. f lT7|) and f JT8|) . Uo(x) and Hi (a;) are given as 

follows: 

T~ fx 2 



no(,H^exp^--), (19) 



= 75 {w - w) exp (-£) + Cexp ' (2o) 

where C is determined by the normalization condition [JdxP^\x) = 1]: 

c 



4V27TQ 3 / 2 
8 



Substituting Eqs. f|T9|) and f )20|) in Eq. f fT6|) . we obtain the stationary distribution Pgp(x). 

Figure |2] shows stationary distributions for the parabolic potential case calculated by 
HFPE (solid lines), FPE (dashed lines), and MC simulation (circles). We also show the RMS 
distance between the HFPE (FPE) and MC distributions in Fig. EJ^a), (b), and (c). We see 
that the densities of the HFPE at x = are higher than those of FPE under the existence of 
noise-intensity fluctuations. Because the stationary distributions with noise-intensity fluc- 
tuations are approximately realized as superposition of many Gaussian distributions with 
different variances, those with smaller variance make the stationary distributions sharper. 
It is important to note that the stationary distributions of the HFPE (solid lines) are not 
Gaussian. This non-Gaussianity is derived from the existence of higher-order derivatives 
than the second in the HFPE (see Sec. 16.11) . From the RMS values, we see that the distance 
between the HFPE and MC distributions is smaller than that between FPE and MC, in- 
dicating a better agreement of HFPE distributions. This result supports the reliability of 
our approximation scheme. The HFPE stationary distribution of 7 = 15 [Fig. EJ^a)] exhibits 
better agreement than the 7 = 5 case [Fig. EJ^c)] because larger 7 yields better approxima- 
tion for the adiabatic elimination technique. The inset in Fig. |2]^a) shows a log-scale plot for 
the x > region, showing that the HFPE stationary distribution has fatter tails, which in 
turn implies that the existence of noise-intensity fluctuations makes distributions fatter (see 
Sec. 16. ip . Furthermore, the HFPE stationary distribution exhibits better agreement with 
the MC simulations than the FPE in tail areas, indicating that our approximation scheme 
offers reliable results even in tail areas for sufficiently large 7 and sufficiently small e. 

We next calculate the stationary distributions of the quartic bistable potential [U(x) = 
x 4 /4 — x 2 /2] driven by additive SIN [g(x) = 1]. Because bistable potentials can represent 
switching dynamics, they are very important in many fields. Hq(x) and H\(x) are given by 



n x (x) = 

with 





FIG. 2: (Color online) Stationary distributions for the parabolic potential with additive noise 
calculated by HFPE (solid lines), FPE (dashed lines), and MC (circles). Parameter values are (a) 
D x = 0.3, D s = 0.2, 7 = 15, and a = 0.5 (Q = 0.135, e = 0.0107); (b) D x = 0.5, D s = 0.3, 7 = 10, 
and a = 0.5 (Q = 0.275, e = 0.0355); and (c) D x = 0.1, D s = 0.1, 7 = 5, and a = 0.1 (Q = 0.011, 
e = 0.00255). The inset in (a) is plotted on a log scale (x > 0) and the others on linear scales. 

where I n (z) is the modified Bessel function of the first kind and C is to be numerically 
evaluated by the normalization condition. 

Figure |3] shows the stationary distributions calculated by the three methods for the quartic 
bistable potential case. The meaning of each symbol is the same as in Fig. [2j We also see 
that the densities of the HFPE at stable sites (x — — 1 and 1) are higher than those of the 
FPE under the existence of noise-intensity fluctuations. The inset in Fig. |3]^a) (a log-scale 
plot) shows that the distribution has fatter tails under noise-intensity fluctuations. As in 
the case of the parabolic potential, the stationary distributions under additive SIN are not 
Boltzmann-Gibbs distributions. According to the RMS values, we also see a good agreement 
between the HFPE and MC results. From the log-scale plot, the HFPE result shows better 
agreement than the FPE also in tail areas for sufficiently large 7 and sufficiently small e. 
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FIG. 3: (Color online) Stationary distributions of the quartic bistable potential with additive noise 
calculated by HFPE (solid lines), FPE (dashed lines), and MC (circles). Parameter values are (a) 
D x = 0.1, D s = 0.1, 7 = 20, and a = 0.4 (Q = 0.026, e = 0.00142); (b) D x = 0.5, D s = 0.1, 
7 = 10, and a = 0.5 (Q = 0.175, e = 0.0157); and (c) D x = 1, D s = 0.1, 7 = 10, and a = 0.1 
(Q = 0.11, £ = 0.0127). The inset in (a) is plotted on a log scale (x > 0) and the others on linear 
scales. 

5.2. Quartic Bistable Potential with Multiplicative Noise 

We consider the quartic bistable potential as in the previous section [U(x) = x 4 /4 — 
x 2 /2], but the system in this section is driven by linear multiplicative noise [g(x) = x]. 
Multiplicative noise plays an important role in many phenomena. Ilo(x) and Ili(x) are 
given by 

no(^ = >l^exp(-|lY (21) 



ni(x) 



C\x\ 



exp 



f_\ _L_ 

2QJ + ZQ 3 



x\ 



exp 



x*\ 
2Qj 



x 



- \{2Q + l)x 4 + ~ (4Q 2 + 6Q + 3) x 2 - log 



x 



11 



where Z and C are 

' 1 
2Q 



Z = (2Q)^r 



C = \ , , f 6 log (2Q) + 6^ - 4Q(Q + 3) - 11 

i2Q3(2Q)wr(^j i \2g; 

Here, ip(x) is the digamma function [ip(x) = d x logT(x)]. We see that n (a;) is composed of 
two Gamma distributions for x > and x < 0. Depending on the effective noise intensity 
Q, n (x) is unimodal (Q > 1) or bimodal (0 < Q < 1). 

Figure 0] represents the probability densities of three typical cases. Figure @Ja) shows 
the bimodal stationary distribution (Q < 1). In this case, the existence of noise- intensity 
fluctuations makes densities around x — 1 and x — — 1 higher. This is also observed in 
the additive noise case. Figure \Mjo) shows qualitatively different shapes between stationary 
distributions of FPE and HFPE. Because Q = 0.96 < 1 in Fig. g(b), the FPE stationary 
distribution (dashed line) exhibits bimodality. On the other hand, the stationary distri- 
bution of HFPE is unimodal (this result agrees with the MC simulation). This indicates 
that the transition from bimodal to unimodal is induced by the existence of noise-intensity 
fluctuations, and the FPE stationary distribution does not correctly reflect this property. 
Figure H](c) also shows the qualitatively different results between FPE and HFPE stationary 
distributions. It is interesting to see that the HFPE stationary distribution is trimodal, 
which is not observed without the existence of noise-intensity fluctuations (the MC sim- 
ulation also exhibits trimodality). This trimodal distribution is a result of superposition 
of unimodal and bimodal distributions. This result shows that the HFPE formulation is 
essential to understand systems driven by multiplicative noise. In all figures, the RMS val- 
ues show that the HFPE distributions provide more reliable results than FPE distributions. 
From the inset in Fig. 11(a) (a log-scale plot), we see that the GPFE distribution agrees with 
the MC distribution for tail areas for sufficiently large 7 and sufficiently small e. 



5.3. Gene Expression Model 

Langevin equations are often used for describing stochastic chemical reactions. In order 
to show that our approximation scheme can be applied to systems including non-trivial drift 
terms, we apply it to a stochastic gene expression model |34j. A simple genetic expression 
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(b) £ = 0.0600 

Q =0.96 

D x =1 
D s = 0. 
7= 20 
a = 0.4 



(a) e= 0.0125 

Q = 0.2 




P st (aO 



HFPE (RMS = 0.00647) 
FPE(RMS = 0.0113) 
MC 




(c) £ = 0.106 

Q =0.925 

= 0.94 
D s = 0.97 
7=9 
a = 0.12 



r P.t(*) 



11FPE(RMS = 0.0143) 
FPE(RMS = 0.0181) 
MC 




FIG. 4: (Color online) Stationary distributions of the quartic bistable potential with linear multi- 
plicative noise [g(x) = x] calculated by HFPE (solid lines), FPE (dashed lines), and MC (circles). 
Parameter values are (a) D x = 0.1, D s = 1, 7 = 20, and a = 1 (Q = 0.2, e = 0.0125); (b) D x = 1, 
D s = 0.8, 7 = 20, and a = 0.4 (Q = 0.96, e = 0.06); and (c) D x = 0.94, D s = 0.97, 7 = 9, and 
a = 0.12 (Q = 0.925, e = 0.106). The inset in (a) is plotted on a log scale (x > 0) and the others 
on linear scales. 



model with self-regulation is given by 



ax 



dx 

dt x 2 + k 



bx + c, 



(22) 



where x is the concentration of a transcription factor (protein), and a, b, c, and k are model 
parameters taking positive values (for details, see Ref. |34j). In this model, the potential 
function is given by 



— -= J H — (a + c)x. 

VkJ 2 



(23) 



The units of a, b, and c are min -1 , implying that the decay time of Eq. (122]) is on the 
order of minutes. Because gene expression goes through various fluctuations, some of the 
fluctuations have faster decay times than x. Then the model can be cast in the form of 
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FIG. 5: (Color online) A potential function U(x) for the gene expression model with a = 6min , 
b = lmin -1 , c = 0.4min , and k = 10. U(x) exhibits bistability with these parameters. 



Langevin equation (j5J with 



/(*) 



(IX~ 



bx + c, 



(24) 



x 2 + k 
g{x) = 1. 

Specifically, we employed a = 6min _ , b = lmin - , c = 0.4 min -1 , and k = 10, as in 
With the adopted parameter values, Eq. (1251) exhibits bistability (Fig. EJ). It is 



Ref. 



considered that this bistability is responsible for the genetic switch. 

Using our approximation method, Hq(x) and Hi(x) are given as follows: 



n (x) 



cxp 



U(x) 



(25) 



ni(x) 



C exp 



U(x) 



8Q 3 Z 



cxp 



U(x) 



-2b 3 x 4 + 8b 2 (a + c)x 3 + 12b (-a 2 - 2ca - c 2 + bQ) 



+8 {a 3 + 3ca 2 + 3 (c 2 - b(bk + Q)) a + c 3 - 36cg} x 
{3xa 2 + 4(bk - 2Q + cx)a - 8Q(c - bx)} 



x 



x 2 + k 



-3aVk (5a 2 + 12ca + 8c 2 - 86 2 A;) arctan {^=j 



+24ab(a + c)Hog (x 2 + k) + 



2afc 



{-kxa 2 + 6A;Qa + 8Q 2 a;) 



(26) 



(x 2 + ky 

where Z and C are normalizing constants, which are numerically evaluated so that 



/ dx n, 



1 and J da;P f 



X 



1. Although the concentration x cannot take nega- 



tive values, the obtained stationary distribution has a very small magnitude in the x < 
areas, which can be ignored in practical calculations. 
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FIG. 6: (Color online) Stationary distributions of the stochastic gene expression model calculated 
by HFPE (solid lines), FPE (dashed lines), and MC (circles). Parameter values are D x = 0.04, 
D s = 1, 7 = 4, and a = 1 (Q = 0.08, e = 0.025). (a) is plotted on a linear scale and (b) is on a log 
scale. 

Figure represents probability distributions calculated by the three methods. The dis- 
tributions in (a) and (b) are plotted on linear and log scales, respectively. We see that the 
density of the HFPE (solid line) at the stable site (x ~ 4.3) is higher than the case without 
noise-intensity fluctuations (dashed line). From the RMS values in Fig. E^a), the HFPE 
stationary distribution agrees well with the MC simulations. In the tail areas [Fig. E^b)], 
the HFPE result exhibits good agreement with the MC simulations, whereas the FPE result 
does not. This result indicates that our approximation method can be applied to systems 
with non-trivial drift terms. 



6. DISCUSSION 



6.1. Effects of higher-order derivative terms 

In the projected time evolution equations of our model, recursive expansion leads to 
equations having an infinite number of terms. In this paper, we have truncated at 0(7 _1 ), 



including derivatives only up to fourth order. Pawula's theorem [35] indicates that to guar- 
antee the positivity of the distribution functions, the truncation must be after the first- or 
second-order term (this corresponds to FPE). Otherwise, the equations must include an infi- 
nite number of terms. However, it has been reported that time evolution equations truncated 
at n > 3 are practically meaningful 36M38| . Indeed our approximation including derivatives 



of higher order than two have shown better results than FPE solutions. In particular, it 
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has been shown in Figs. E^b) and (c) that stationary distributions of FPE and HFPE have 
qualitatively different shapes, and our results using HFPE are supported by MC simulations. 
Ref. 39] gives an intuitive explanation of the effect of higher-order derivatives, which are 



not included in the conventional FPE, on the positivity of solutions. Following Ref. 
we explain the parabolic potential case. Using Eq. (TT3"j) . the FPE of the parabolic potential 
case is given by 



i^Hs'+O^- (27) 



Let x m be points where P(x; t) is locally minimal with respect to x. According to require- 
ments of the minima and the positivity of P(x m ; t), d t P(x m ; t) is always positive: 



d d 
—P{x m -t) = x—P{x-t) 



2 



+ P(x m ;t) + Q—P(x;t) 



> 0. (28) 



Eq. ( f28j) guarantees the positivity of initially positive solutions because at the minima 
P(x m ; t) is increasing as a function of t. On the other hand, the positivity of dtP(x m ; t) is not 
generally satisfied in HFPE cases, since the fourth-order derivative term (R/~/)d%.P(x; t)\ x=Xm 
can take any values. However, the effect of the fourth-order term on the sign of d t P(x m ; t) 
is negligible under sufficiently smooth initial distributions and sufficiently large 7, since the 
fourth-order term is on the order of r y~ l . From the above explanation, it is generally expected 
that solutions of HFPE can satisfy positivity under sufficiently large 7. 

In order to analyze effects of higher-order derivatives on moments for the additive noise 
case, we consider the following moment expansion using Eq. ( 1121) : 



^- ( x n ) = n (f(x)x n - 1 ) + Qn(n - 1) <x n ~ 2 ) + -n(n - l)(n - 2)(n - 3) (x n ~ 4 ) , (29) 
at 7 

where (• • ■ ) represents the expectation with respect to P(x;t) [(A(x)) = f dx A(x)P(x; t)]. 

We consider the stationary case of the parabolic potential case [f(x) = —x\. Equation ( 1291) 

in this case is described by 

<ar> = 0, (x 2 )=Q, (a; 3 ) = 0, <x 4 > = 3Q 2 + ^. 

Kurtosis is given by 

6D.«c? + D) 

7—>oo yields the Gaussian stationary distribution with k = 0. Equation (I30p shows that 
smaller 7 yields a stationary distribution with larger kurtosis. This means that the distri- 
butions have fatter tails under noise-intensity fluctuations. This result agrees with those of 
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Sec. 15.11 In Ref. |13|], fat-tailed stationary distributions are derived not from time evolution 
equations but from a calculation of the expectation with respect to a noise-intensity distri- 
bution. Our results show that derivatives of orders higher than two are required to yield 
fat-tailed stationary distributions in the parabolic potential. Furthermore, the importance 
of the fourth-order derivative was pointed out in Ref. [20|, which considered superstatis- 
tical Brownian particles using mesoscopic nonequilibrium thermodynamics 40[. Ref. 20 1 
obtained a time evolution equation in the Fourier space. Inclusion of a quartic term in the 
terms of a Fourier conjugate variable indicates the existence of a fourth-order derivative 
in real space. In specific cases, stationary distributions of their model are very close to a 
g-Gaussian distribution, which are strongly connected to Beck's superstatistics. 

The stationary distributions of the FPE (Sec. 15. ip belong to the exponential family, which 
are strongly connected to the Boltzmann-Gibbs-Shannon (BGS) entropy. On the other 
hand, the stationary distributions of the HFPE with fatter tails are not exponential. Non- 
exponential distributions in statistical mechanics have been extensively discussed in Tsallis 



statistics 
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Il6| . It was shown that long-range interactions or long-time memory effects are 
responsible for such distributions. Since the noise intensity of our model is governed by the 
Ornstein-Uhlenbeck process, correlation time of the noise intensity serves as the memory 
effect, which is not the case for the white noise. 



6.2. Relation to colored noise 



As argued in Sec. HJ colored noise is a generalization of white noise and also incorporates 
mesoscopic time-scale inhomogeneity. However, many approximation schemes for colored 
noise, including adiabatic elimination, can take into account a time-correlation effect [i.e., 
0(t)] within the conventional FPE. The necessity of derivatives of higher order than two 
in our systems indicates that systems driven by SIN are essentially different than those 
governed by colored noise. 

For colored noise driven systems, one popular approach to obtain approximate FPE is the 
unified colored- noise approximation (UCNA), which is based on the adiabatic elimination 
approach 26]. In Ref s. 4lLl42j]. UCNA was confirmed as a reliable Markovian approximation 



by the path-integral. In Refs. 41 



42], they derived FPE using a Markovian Lagrangian 



function in the path-integral representation, where non-Markovian terms, such as x and 



17 



FIG. 7: (Color online) Stationary distributions of the quartic bistable potential with (a) additive 
and (b) multiplicative noise for 7 = 0.1 (squares), 1 (circles) and 10 (triangles). Parameter values 
are (a) D x = 0.5, D s = 0.1, and a = 0.5 (Q = 0.175) and (b) D x = 0.94, D s = 0.97, and a = 0.12 
(Q = 0.925). These data are from MC simulations with iV = 5 x 10 7 samples. Solid lines are 
included ELS cL guide to the eye only. 

x n (n > 2), were ignored. The path- integral was applied to non-Gaussian colored noise in 
Ref. (43)], and an effective Markovian FPE was obtained in this case. Since our obtained 
equation includes derivatives of higher order than two, a naive application of the path- 
integral techniques used in the colored noise case may not be possible. However, it is 
important to understand our approximation scheme from the viewpoint of the path-integral, 
because the path-integral has an intuitive interpretation in stochastic processes. We leave 
these for future studies. 

6.3. Cases of the other 7 range 

In preceding sections, we described the approximation scheme for the 7 ^> 1 case. It 
is important to consider the other 7 range. In this section, we describe MC simulations 
carried out for 7 = 0.1, 1, and 10 for two settings: the quartic bistable potential with 
additive and with multiplicative noise. Parameters (D x , D s , and a) of the additive and 
multiplicative noise cases are identical to those shown in Figs. EJ^b) andHJ^c), respectively. 
These results are shown in Fig. [7J Fig. [7(a) shows that higher peaks emerge at x = ±1 for 
smaller 7. For Fig. [7(b), densities at x — —1, 0, and 1 are higher for smaller 7. Similar 
trends have also been observed in Sees. 15 . 1 1 and [5T2l Results of Sec. 15.21 and this section show 
that trimodality is induced when the noise-intensity fluctuations are meso/macroscopic and 
the fluctuation width is large. The stationary distributions of smaller 7 can be explained 
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TABLE I: RMS distance between HFPE (or FPE) and MC distributions for parabolic and quartic 
bistable potentials with additive noise, quartic potential with multiplicative noise, and a gene 
expression model. 





parabolic potential 
with additive noise 


quartic bistable potential 
with additive noise 


Fig. 2(a) 


Fig. 2(b) 


Fig. 2(c) 


Fig. 3(a) 


Fig. 3(b) 


Fig. 3(c) 


HPFE 


0.00418 


0.00788 


0.0934 


0.0120 


0.0173 


0.0253 


FPE 


0.0184 


0.0229 


0.152 


0.0338 


0.0267 


0.0329 





quartic bistable potential 
with multiplicative noise 


gene expression 
model 


Fig. 4(a) 


Fig. 4(b) 


Fig. 4(c) 


Fig. 6 


HPFE 


0.00683 


0.00647 


0.0143 


0.00446 


FPE 


0.0165 


0.0113 


0.0181 


0.0163 



by the superposition of different distributions. For 7 — > 0, the case reduces to that of 
superstatistics, where the stationary distribution is given in a Bayesian fashion: 



PJx) 



dsPjx\s)P(s). 



Here, P st (x\s) is the stationary distribution taking s as a parameter, and P{s) is the distri- 
bution of s. 

Ref. [18| considered a different superstatistical model of multiplicative noise processes: 



which assumes that multiplicative noise exponent a fluctuates [£(£) is white noise]. Station- 
ary distributions of their model also exhibit the transitions of stationary distributions. 



7. CONCLUDING REMARKS 

In this paper, we have derived the time evolution equation of systems driven by SIN 
using adiabatic elimination under the assumption 7 ^> 1. The obtained HFPE is of 0(7 _1 ) 
and contains derivatives up to fourth order. We have calculated the stationary distributions 
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of the obtained equation by using the perturbation expansion and applied them to three 
different cases. Table [J summarizes the RMS distance between the stationary distributions 
calculated by the HFPE (or FPE) and MC which have been shown in Figs. (2J El Eland El It 
clearly shows that the results of the HFPE are better than those of FPE. We have pointed 
out that SIN makes densities at stable sites higher. Furthermore, SIN induces the transition 
from unimodal to bimodal in the quartic bistable potential driven by multiplicative noise. 
For specific parameters, the stationary distribution exhibits three peaks, which are not 
observed in conventional white noise cases. Because the stationary distributions obtained 
by Eqs. (fT7|) and f lTS]) are given by the first-order differential equations, their analytical 
expressions can be calculated for many general configurations. The applicability has been 
shown with the stochastic gene expression model, which has a non-trivial drift term. This 
indicates that our method can be applied to many real-world phenomena. Furthermore, 



higher-order derivatives in HFPE have been given much attention in recent years [39|, |44 |. 
Analysis of higher-order terms in our systems are also important, and so we intend to pursue 
this in future studies. 
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Appendix A: Derivation of the Higher-Order Fokker Planck Equation 

adiabatic elimination with eigenfunction 



In this appendix, we explain the procedures o: 
expansion used in our model, following Refs. 23|, 



33]. 



We first expand P(x, s; t) into the complete set ip n (s) of the operator L s (s) [Eq. fl9])]: 

oo 

P(X, S-,t) = J2 P n& t)Ms), (Al) 
n=0 

where ip n {s) are eigenf unctions of the following equation: 

L a (s)V> n (s) = -\ n il) n (s). 
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Eigenvalues X n and eigenfunctions if) n (s) are given by 

A„ = n, (A2) 



Ms) = \I^d s ^\ H M^ (V) > (A3) 



where 77 = \Jl/ (2D s )(s — a), and H n (rj) is the nth Hermite polynomial. We introduce the 
adjoint function Vn( s ) °f V'n(s): 

^ = ^(77), Vj(*) = l- (A4) 
The ort honor mality and complete relations are 

(ViWmOO) = S nm , (A5) 

where (• ■ ■ ) denotes an integration with respect to s [(A(s)) = J ds A(s) ). Using Eqs. ()A4j) 
and 0A50 . P(x,t) (= (P(x,s,t)) ) is given by 



P(x;t) = ^P„(z;*) (^J(s)V'n(s)) = ^ofo*)- (A6) 

n=0 

By multiplying both sides of Eq. (J6]) by Vm( s ) an d integrating out s, we obtain [P m = 

P m (x;t)] 

fs a 00 

di Pm = ~dx~ fPm + E>X ^ i^m^n) A g P n - 7 A m P m , (A7) 

n=0 

where A g is an operator defined by Eq. (ITT1) . 

Using relations of the Hermite polynomial, specifically, if n+ i(,2) = 2zH n (z) — 2nH n -i(z) 
and Hq(z) = 1, the following relation holds: 

(ipl l (s)s 2 '4) n (s)) = 2D s m(m - l)5 m - 2 ,n + 2^/2Z> s oim5 m _i in + -D s 5 m+2 , n 

+y/2D~ a a6 m+ljn + |2D S fm + + a 2 j 5 m , n . (A8) 
According to Eqs. (jA7B and (]A8[) . time evolution of Po is governed by 

dt P ° = ~dx~ fP ° + + ^ AgP ° + >fiD> aD * A 9 P i + ^D s D x A g P 2 . (A9) 



Equation AA9D includes Pi and P 2 , although we want to obtain a closed equation of P . We 
calculate P m (m > 1) by Eq. (1X7]) : 

Pm = irrfl"^ + ^E(W A » P 4 (form>l). (A10) 

Am7 I °^ n=0 J 
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In 

up 



Eq. (1A10j) . we ignore a time derivative term. By substituting Eq. ( lAlOj) into Eq. (1A9j) . 
to 0(7" 1 ), we obtain Eq. §W}). 
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